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1. Introduction 

The production of hadrons in e + e~ annihilation is one of the best-studied processes in high- 
energy physics. It played a crucial role in establishing the correctness of QCD as the theory 
of strong interactions. Studies of hadronic final states in e + e~ annihilation at PEP, KEK, 
PETRA, SLD and LEP were instrumental for understanding jets, advancing perturbative 
QCD computations, developing parton showers and investigating non-perturbative QCD 
effects in high-energy collisions. 

At LEP experiments, exclusive processes with up to five jets in the final state were 
studied in detail; inclusive measurements are available for up to six jets [1—16]. Theoreti- 
cally, e + e~ annihilation cross sections into two [17, 18] and three [19-21] jets are known at 
next-to-next-to-leading order (NNLO) in perturbative QCD, the production of four jets is 
known at next-to-leading order (NLO) [22-26], while only leading order (LO) predictions 
were available so far for five-jet observables. Resummed results for jet rates have been also 
obtained [27-31]. 

High-quality data, primarily from LEP measurements at the Z-pole, open up an op- 
portunity to perform very accurate studies in jet physics [32]. Since, by now, gross features 
of QCD are well understood, the interest shifts towards subtle details which can be revealed 
only through dealing with complicated final states and improving the accuracy of theoret- 
ical predictions. Because the n-jet rate is proportional to the strong coupling constant at 
high power <r n j e t ~ a s -2 j leading-order predictions for n-jet observables for n>2 are very 
uncertain and an improved theoretical description of such final states is desirable. When 
improved descriptions become available, they are used to study interesting properties of 
hadronic final states including event rates, event shapes and ultimately to extract the value 
of the strong coupling constant a s [18-20,23,33-38]. 
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Infra-red safe observables, traditionally studied in e + e~ annihilation, are dominated 
by short-distance physics; they are computable in perturbative QCD up to corrections 
suppressed by inverse powers of some large energy scale related to e.g. the average relative 
transverse momentum of jets in a given process. Those power (or hadronization) corrections 
are typically estimated using event generators such as PYTHIA [39], HERWIG [40], and 
ARIADNE [41], under the assumption of a complete factorization of non-perturbative and 
perturbative physics. This factorization implies that the hadronization corrections to an 
infrared-safe observable O can be estimated as 

H l [0] = 5b*-, (1.1) 

^part 

where C^adr anc ^ O part are the values of the observable O computed at the hadron and at 
the parton level with the event generator i. Because hadronization corrections are assumed 
to be factorizable, they can be used to "improve" the perturbative prediction pt for the 
observable under study. Hence, one defines the quantity Oi mpr = H l [0] O pt , and compares 
it to experimental data 1 . 

While the use of the procedure that we just described is widespread, it is clear that 
it can not be fully valid. Indeed, one can imagine that, for a particular observable O, the 
event generator i happens to reproduce its measured value, O^adr = OData- As a result, 
^-^impr reads 

Q ^ 

Oimpr = Ooata • (1-2) 

^part 

Clearly, Oi mpr can only be equal to Od ata if Opt = Opart > but this equality can not hold true 
for a variety of reasons, including different approximations in parton/dipole showers and 
fixed-order perturbative calculations, different dependencies of O pt and O part on the renor- 
malization scale fi, etc. Differences at the perturbative level are particularly worrisome 
when considering high-multiplicity final states, since standard event generators routinely 
used in e + e~ studies are based solely on e + e~ —¥ qq and e + e~ — > qqg matrix elements. 
This implies that, since event generators are tuned to data, hadronization corrections as 
defined by eq. ( |0| ) contain both non-perturbative and perturbative effects; the latter are 
present to compensate for deficiencies of the partonic part of a particular event genera- 
tor 2 . To the extent that both perturbative and hadronization corrections are small, the 
inconsistency of the whole procedure may not be very apparent, but it becomes evident 
if those corrections are large. It is very likely that these issues are important for five-jet 
production at LEP. Indeed, because e + e~ — > 5 jets involves a high power of a s , NLO QCD 
corrections are expected to be large. In addition, a correct description of five hard, well- 
separated partons is difficult for traditional event generators, so that all the problems of a 
conventional approach to estimating hadronization corrections can be exposed by studying 
five-jet observables. 

1 There are alternative approaches to estimate hadronization corrections that address both theoretical [42] 
and experimental [16,43] aspects of this procedure. 
2 See ref. [44] for a related discussion. 
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By extracting hadronization corrections with traditional event generators - PYTHIA, 
HERWIG and ARIADNE - we find these corrections to be large and generator-dependent. 
This is unfortunate since it implies a large spread of "improved" predictions for five-jet 
observables when perturbative and hadronization effects are combined. For this reason, 
we believe it is important to obtain hadronization corrections from an event generator 
whose perturbative part is up to the task of describing the production of five hard, well- 
separated partons. Within the context of event generation, the implementation of such a 
description requires a consistent matching between high-multiplicity matrix elements, and 
parton/dipole showers. In this paper, we employ the SHERPA event generator [45,46], 
which implements the CKKW [47] matching prescription 3 . By calculating hadronization 
corrections with SHERPA/CKKW, we find that they are relatively small (see the right pane 
of fig. ^), in particular in the kinematic region where perturbative QCD is reliable. This 
is what we expect since, with CKKW matching, the perturbative description of five-jet 
production provides a good approximation to the actual physical process. As a conse- 
quence, when we use SHERPA/CKKW to extract non-perturbative corrections according 
to eq. the results are less contaminated by perturbative contributions, compared to 

the case when traditional event generators are used for this purpose. Since these consid- 
erations apply to hadron collisions as well, our results have obvious implications for jet 
physics at the LHC. 

The goal of this paper is to investigate five-jet production at LEP, using theoretical 
predictions that are accurate at NLO in QCD. We discuss hadronization corrections and 
show that they depend significantly upon the event generator that is used to estimate 
them. We extract the strong coupling constant by fitting the NLO QCD predictions for 
the distributions of the five-jet resolution parameter 2/45 and the five-jet rate R$. The 
remainder of the paper is organized as follows. In Section |2| we discuss the technical details 
pertinent to the computation of the NLO QCD corrections to e + e~ — > 5 jets. In Section || 
a phenomenological analysis of five-jet production at LEP is reported. In Section |] the 
value of the strong coupling constant is extracted. In Section [| we present our conclusions. 
Details of the fit procedure that we use in our analysis of the strong coupling constant are 
described in the Appendix. 

2. Technical details 

The computation of NLO QCD corrections to any process or observable in the context of 
a subtraction formalism involves the evaluation of the following contributions: the one- 
loop virtual corrections, the real-emission corrections and their subtraction terms, and the 
finite remainders obtained from the analytical integration of the subtraction terms over 
the degrees of freedom of the unresolved parton(s). The one- loop virtual amplitudes for 
e + e~ — > 5 jet required in this paper are computed within the generalized D-dimensional 

3 We emphasize that it is not the parton/dipole shower per se that makes the difference, but the possibility 
to match it to high-multiplicity matrix elements. We expect to obtain quantitatively similar results with 
HERWIG, PYTHIA, and ARIADNE if these event generators are supplemented with a matrix-element- 
matching procedure. 
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unitarity framework, as described in Refs. [48-50]; a few technicalities relevant to the five-jet 
case are given in sect. ^4]. The remaining contributions are calculated using MadFKS [51], 
which is also employed to perform the integration over the phase space of the final short- 
distance cross sections. Details of MadFKS relevant for this computation are reviewed in 
sect. 



2.2 



2.1 Calculation of one-loop amplitudes 

Within the context of generalized D-dimensional unitarity, we compute the so-called prim- 
itive amplitudes [52], which are gauge-invariant subsets of color-ordered amplitudes. The 
color decomposition of amplitudes that we need in this paper coincides with the color de- 
composition of QCD amplitudes without any colorless vector boson ~f*/Z. We use the color 
decomposition introduced in ref. [53]. The relation between primitive and color-ordered 
amplitudes that we need in this paper can be found in ref. [50]. 

In general, the amplitudes needed for the NLO calculation of e + e~ — > "y*/Z — > 5 
partons are very similar to the amplitudes with a VF-boson and five partons. The main 
differences between amplitudes that involve charged and neutral currents originate from 
the different couplings of W^-bosons and 7* jZ bosons to fermions. In particular, in the 
Standard Model, the coupling of the W-boson to fermions is only left-handed, while for 
*y*/Z one has to sum over left- and right-handed states. Indeed, the Z-boson couples to 
fermions via a vector and a vector-axial coupling gvj^" + 9A7^75- We can rewrite this 
coupling through left- and right-handed projection operators -Pl,r = (1 T 7s)/2 as 

9vY + QA^ls = <?L7 M ^L + QrI^Pk, (2.1) 

where gi,R = 9V-F9A- Therefore, all the one-loop amplitudes relevant for this paper can be 
obtained by considering quarks coupled to a vector current only; the only subtlety is that 
the vector current should couple to left- and right-handed quarks with different strength. 

As we already mentioned, many of the non-trivial amplitudes that we need are iden- 
tical to the amplitudes calculated for — > W + 5 partons case [50]. We note, however, 
that new amplitudes appear if 7* jZ couples directly to a loop of virtual fermions. There 
are two reasons for this. First, such amplitudes are not present in the W-boson case stud- 
ied in ref. [50] because of charge (or flavor) conservation. Second, certain parts of those 
amplitudes are related to an axial anomaly and, therefore, violate the symmetry between 
the vector current and the axial- vector current. It is interesting to remark that within the 
context of generalized D-dimensional unitarity, a correct computation of the axial anomaly 
entails a literal implementation of the 't Hooft-Veltman prescription for dealing with 75 in 
dimensional regularization. We have calculated the required amplitudes where j*/Z cou- 
ples to the fermion loop but we neglect them in this paper since previous experience with 
those amplitudes shows that they are very small [22,23,54,55], especially when compared 
to the residual theoretical uncertainty of the five-jet rate. 

The cross section involves Born and virtual amplitudes with two or four quarks in the 
final state. We relate amplitudes involving a Z-boson exchange to amplitudes that only 
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involve the photon (vector) exchange: 



A z (q\,q-\,e+e_ „,.••) = Qn Qe ^z(s) 



9R9R^RxSR p ,A 1 *(q_,q + ,e-,e + , . . .} 



+ 9R9l^R\^LtiA r (q^,q + ,e + ,e-, . . .) + g q L g R 8 L \5 RlJ ,A 7 *(q + , g_, e_, e + , . . .) 



+ 9 l 9l^l\8l^* (q + ,q..,e + ,e-,...) 



(2.2) 



In eq. ( |2.2| ) all particles are outgoing, the dots denote all gluon momenta, and the propa- 
gator factor V(s) is given by 

s - M| + iTzMz 

where Mz and are the mass and the width of the Z boson respectively. The left- and 
right-handed coupling of the electrons (/ = e) and quarks (/ = q) to the Z boson are 
explicitly given by 

J_J f _ 2(Tj - Qf sin2e w ) f _ f f _ -2Qf sm26 w 

9l-9 v -9a- —(20^) > 9 R -9v + 9 A - sin(2 e w ) ' {2A > 

where 6 W is the Weinberg angle, T 3 e = T 3 d = -1/2, Tg = 1/2 are the values of weak isospin 
for quarks and leptons and Q e = —1, Q u = 2/3, Q d = —1/3 are the respective electric 
charges. 

The amplitudes involving four quarks and a Z-boson can be written in a similar way. 
The amplitude that involves different quark flavors can always be written as the sum of 
two amplitudes, where the Z-boson couples to a particular quark line 

Az(q\,q-x, Qa, Q-a, e+, eI M , . . .) = A z (q\, q~\, Qa, Q-a, rf, • • •) 

+ A z {Qa, Q-a, q\, q~x, e+,eZ^, . . .)■ ( 2 -5) 

Note that when amplitudes A are computed, the Z-boson is allowed to couple to the quark 
flavor indicated by first and second argument of A. The expression for Az amplitude in 
terms of A* is identical to eq. (|2.2|), so we do not repeat it here. Finally, if the flavor of 
the two quark lines is the same, we include the symmetry factor 1/4 and anti-symmetrise 
with respect to the exchange of the quarks or anti-quarks. We point out that the full 
one-loop matrix elements squared that we use in this paper were checked against a similar 
computation performed by the BlackHat collaboration [56-58], and complete agreement 
was found. 

2.2 Real emission corrections with MadFKS 

MadFKS is based on the (FKS) subtraction formalism of ref. [59]. The implementation of 
the FKS procedure is fully automated in MadFKS. In essence, MadFKS goes through the 
following steps. First, it determines the partonic processes that contribute to a given phys- 
ical reaction, and their singularity structures. Then, it constructs the real-emission matrix 
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Process 



# of FKS pairs 



e + e- -> qqgggg 3 

e+e~ -)■ qqq'q'gg 7 

e + e~ -> g<M#c/ 4 

e + e~ — > qqq'q'q"q" 3 

e + e~ — )• qqqqq'q' 2 

e + e~ — )■ qqqqqq 1 



Table 1: Numbers of FKS pairs for the various real-emission processes that contribute to the 
five-jet cross section. See the text for details. 



elements, their subtraction terms, the finite remainders, and the Born matrix elements. Fi- 
nally, it proceeds to the actual computation, by sampling (possibly with multi-channeling 
techniques) the phase space, by evaluating the short-distance cross sections, and by return- 
ing weighted parton-level kinematic configurations, that can be used to construct as many 
observables as one likes. All matrix elements, except those of virtual origin, are obtained 
by calls to MadGraph [60] routines. The virtual matrix elements are on the other hand 
computed as described in sect. 2.1. We point out that MadFKS gives, for each phase-space 
point, a four- momentum configuration as input to the code of sect. 2T, which returns three 
numbers - corresponding to the double and single (IR) pole residues, and to the finite part; 
the talk-to between the two codes uses the Binoth-Les Houches interface [61]. As clarified 
in ref. [51], MadFKS integrates quantities that are locally finite in the phase space, and 
in four dimensions. Therefore, the pole residues provided by the virtual amplitudes are 
used only to check that they match those that are known analytically from the subtraction 
procedure, in this way ensuring that KLN cancellation does indeed take place. The only 
output of the code of sect. 2A used in the integration of the short-distance cross sections 
is thus the finite part, defined according to the conventions given in the Appendix B of 
ref. [51]. 

We have performed the calculation of the five-jet cross section in a straightforward 
manner. We included all partonic processes, and performed explicit sums over colors and 
helicities, except in the case of virtual amplitudes. For the latter, we have computed 
separately the leading- and subleading-color contributions, and performed the sum over 
helicities using Monte Carlo methods. We have used five massless quark flavors. We 
remind the reader that the FKS formalism is particularly efficient in keeping the number 
of subtraction terms to a minimum. Furthermore, the real-emission matrix element minus 
the subtraction terms is re-organized into a sum that gets as many contributions as the 
subtraction terms themselves, and that are separately finite, which implies that they can 
be (and are) integrated independently from each other. Physically, these contributions 
corresponds to pairs of particles (called FKS pairs) that can give one soft and/or one 
collinear singularity at most. We report in table [l] the number of FKS pairs relevant to 
the various real-emission processes that contribute to the five-jet cross section. 
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With five massless flavors, the number of independent partonic subprocesses that con- 
tribute to the five-jet cross section is 25. Using the entries of table |l[ this implies 81 FKS 
pairs in total, i.e. 81 independent integrations. On the other hand, the complexity of the 
kinematics is such that even in the context of an adaptive integration it may be very dif- 
ficult to map correctly all the peaks of the Feynman diagrams, and thus to have a stable 
numerical behavior. We have therefore preferred to adopt a multi-channelling integration 
strategy, that in MadFKS follows the same procedure as in MadGraph [62]. In doing so, 
the numbers of integration channels we deal with at the real-emission and virtual level are 
equal to 3620 and 2 x 1408 respectively (the factor of two in the virtual amplitudes being due 
to the independent integration of the leading- and subleading-color contributions). These 
numbers are much larger than the 81 FKS pairs we started with; however, the Feynman- 
diagrams peaks can now be mapped accurately by the integration routines, and relatively 
small statistics is sufficient in each channel to obtain numerical stability. We conclude 
by stressing that the MadFKS integration channels are fully independent. Furthermore, 
they are not determined dynamically (e.g. by performing a preliminary integration of the 
cross section), but are defined a priori, by considering the topologies of the Feynman di- 
agrams that contribute to the relevant partonic processes. The whole organization of the 
calculation is therefore inherently parallel. 

3. Phenomenology of five-jet production 

In this Section we present the results of our calculation. We consider e + e~ — > jets and define 
jets using the Durham jet algorithm [28] with resolution parameter y cu t- The following 
distance between each pair of particles is used in the Durham jet algorithm 

Vij = — (1 - cos dij) , 

where s is the center-of-mass energy of the collision squared, Ei is the energy of the parton 
i, and 9{j is the relative angle between the partons i and j, in the e + e~ center-of-mass 
reference frame. The pair of particles with the smallest distance is clustered together by 
adding their four- momenta, as long as yij < y C ut, and the procedure is then iterated. When 
all distances yij are larger than y C ut> the recombination stops and the number of jets in the 
event is defined to be equal to the number of (pseudo)-particles left at that stage. 

In this paper we consider two observables which we define with the Durham jet al- 
gorithm. The first observable is the differential distribution with respect to the five-jet 
resolution parameter 2/45, normalized to the total cross section for e + e~ —> hadrons, <7tot 
(which we compute at the NLO, i.e. at 0(a s )). The resolution parameter 2/45 is the maxi- 
mal value of y cn t such that a given event is classified as a five-jet event by the Durham jet 
algorithm. We note that 



(3.1) 



where o"f nc j et is the inclusive five-jet production cross section in e + e annihilation. The 
second observable that we study is the five-jet rate RniVcat)- It is defined as follows 

5— jet/ \ 

R 5 (y m t) = acxd [Vcut) , (3.3) 

Ctot 

where (y C ut) is the exclusive five-jet production cross section. It is calculated by 

applying the Durham jet algorithm to the given event, and by requiring that exactly five 
jets are reconstructed. 

When we compute a^da/dlny^ and R§ in perturbative QCD, we obtain a power 
series in the strong coupling constant 

a ' ol d^ = (^r) 3 Ai5im5) + (^r) 4 (^5(2/45) + 36 A 45 (y 45 ) In ^) ,(3.4) 



flsGfcut) = (^-gpj As(ycut) + ^^rj ^ 5 (y cu t) + 36o^ 5 (ycut)ln^J , (3.5) 

where is the renormalization scale, bo = (33 — 2nj)/3 and ray = 5 is the number of quark 
flavors that we treat as massless. The top quark is considered to be infinitely heavy and is 
completely neglected in our computation. It is important to emphasize that the coefficients 
^45,5 and -645,5 depend on j/45 and y cu t, respectively, but not, say, on the total center-of-mass 
energy squared. This feature is a consequence of the following approximations employed 
in our computation: 1) all particles, except the Z boson, are treated as massless; 2) the 
observables that we are interested in are sufficiently inclusive so that the vector and the 
axial currents do not interfere; 3) we neglect triangle fermion diagrams that lead to the 
axial anomaly so that (for equal couplings) vector and axial current contributions to the 
final result are equal. These three points are sufficient to ensure that ^45,5 and .645,5 are 
independent of the electroweak parameters and the center-of-mass energy squared. 

Experimentally [16], five-jet observables are computed using the reconstructed mo- 
menta and energies of charged and neutral particles. Measurements are corrected for 
detector effects, so that final distributions correspond to stable hadrons and leptons, and 
for initial- and final-state photon radiation, which is a sizable correction for LEP2 data. 
Above the Z peak, relevant backgrounds are subtracted; the most important among them 
is PF-pair production. The experimental uncertainties are estimated by varying event- and 
particle-selection cuts. They are below 1% at LEP1 and slightly larger at LEP2. Further 
details of the experimental analysis can be found in ref. [16]. 

In fig. |], we compare ALEPH LEP1 data [16] for 1/crtotdcr/dlny^ 1 with the hadron- 
level predictions of three event generators - PYTHIA, HERWIG and ARIADNE. We 
observe that these event generators describe experimental data fairly well; differences 
between data and theoretical predictions are below twenty five percent in the central 
4.5 < lny^ 1 < 9 region of the distribution, where the statistical accuracy of the data 
is good. This is an impressive accomplishment since a^dcr/dlny^ 1 changes by three or- 
ders of magnitude in this range of In . On the other hand, it is clear from the upper pane 
of fig. H that hadronization corrections are very large and change from 0.5 to 1.5 in that 
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Figure 1: ALEPH data [16] for the y 45 distribution at LEP1, compared to PYTHIA, HERWIG 
and ARIADNE results. The upper panes show detector and hadronization corrections, respectively. 
The lowest pane shows the relative difference between data and event generator predictions. This 
figure was provided to us by H. Stenzel. 



range of lny^ 1 . In addition, it follows from fig. [l] that the difference between hadronization 
corrections, as calculated using different event generators, can be as large as 20-30%. 

We attribute these features to the inability of PYTHIA, HERWIG and ARIADNE 
to describe hard perturbative radiation correctly. Indeed, these programs generate high- 
multiplicity final states starting from hard low-multiplicity processes; they produce addi- 
tional jets by means of parton/dipole showers. Since these showers describe hard large-angle 
emissions only approximately, the so-called hadronization corrections attempt to correct 
for this (perturbative) deficiency. While this problem is unavoidable if traditional event 
generators are used to describe high-multiplicity final states, techniques exist to match 
parton showers and high-multiplicity matrix elements in a consistent manner, thereby im- 
proving the pure-perturbative part of event generators. One such technique is the CKKW 
matching procedure [47], which is implemented as default in the SHERPA event gener- 
ator. The comparison of ALEPH LEP1 data with SHERPA predictions, as well as the 
hadronization corrections derived from SHERPA, are shown in fig. ||. Two hadronization 
models - Lund string [63] and cluster [64] - are employed. In the central part of the 
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Figure 2: ALEPH data for the 2/45 distribution at LEP1, compared to SHERPA results. Two 
hadronization models - Lund string [63] and cluster [64] - are employed. The lower pane in the left 
plot shows the relative difference between Sherpa predictions with the two hadronization models, 
and ALEPH data. In the right plot, the hadronization corrections for the two models are shown. 

distribution, SHERPA results agree with ALEPH data to 20 — 25%, similar to traditional 
event generators. Moreover, in the region of moderately small values of lny^g , where fixed- 
order perturbative description is reliable, the hadronization corrections are below twenty 
percent, in sharp contrast with estimates of hadronization corrections based on PYTHIA, 
HERWIG and ARIADNE. It is important to emphasize that, although in that region of 
In yjg 1 traditional event generators provide slightly better description of data compared to 
SHERPA, this does not mean that hadronization corrections extracted with the former 
codes are more reliable. Indeed, traditional event generators achieve agreement with data 
at the price of very large hadronization corrections. This feature precludes a clear sep- 
aration between long- and short-distance phenomena, which is crucial for the procedure 
outlined below eq. ( |1.1| ) to be meaningful. 

The ALEPH data exhibit a characteristic turnover shape. This turnover means that 
for small values of 2/45, the result is dominated by exclusive five-jet production with very 
small resolution parameter, where fixed order perturbation theory fails and a resummation 
is required to achieve meaningful results. A resummation of a"L 2n and a™L 2n_1 terms, 
where L = In y~^ t , was performed for R§ in [28] , while no resummation is currently available 
for the five-jet resolution parameter distribution. However, there seems to be no region in 
L where this resummation can be valid since two conditions L> 1 and a s L <C 1 should 
be satisfied simultaneously. Taking a s ~ 0.15 as a typical value of the strong coupling 
constant 4 , we find that L should be smaller than 6. On the other hand, practical experience 

4 We take 5 — 20 GeV as a reasonable estimate of the scale of the strong coupling constant for five-jet 
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Figure 3: ALEPH LEP1 data compared to leading and next-to-leading order predictions in QCD, 
without hadronization corrections. We use a s (Mz) = 0.130 at the leading and a s (Mz) = 0.118 at 
the next-to-leading order in perturbative QCD. The renormalization scale is chosen to be O.SMz- 
The uncertainty bands are obtained by considering the scale variation 0.15 Mz < [i < 0.6 Mz- 
Solid lines refer to NLO QCD results evaluated with n = 0.3Mz- 



with resummations suggests that L ^> 5 is what can be considered as a large logarithm. 
Clearly, 5 <C L < 6 leaves very little room for the validity of this approach. It should be 
possible to improve on the resummation by including sub-leading logarithms and matching 
to NLO QCD computations. However, since we do not perform any resummation in this 
paper, we require lny± 5 ,lny~ ut ^ 6 for the comparison of the NLO QCD computation with 
data. Interestingly, a similar upper bound on In y~^ appears because we neglect the mass of 
6-quarks in our computation. This implies that the resolution parameter times the center 
of mass energy should be larger than the 6-quark mass, i.e. sy^ > ml, which translates 
into ln(y^) < ln(s/rii|) £ 6, for s = M§. 

When fixed-order perturbative QCD calculations are compared to experimental data, 
the choice of the renormalization scale becomes an important issue. Traditionally, multi- 
jet observables in e + e~ annihilations are computed in perturbative QCD by evaluating 
the strong coupling constant at the center-of-mass energy. However, for large numbers of 
jets this choice should be reconsidered, since the hardness of each jet decreases with their 
number. Dynamical renormalization scales used in event generators account for this effect 
by relating the choice of the renormalization scale to the event kinematics. Our choice 
of the renormalization scale is also motivated by dynamical considerations. To this end, 
we consider the clustering history of five- and six-parton configurations that results from 
using the Durham jet algorithm. We compute the average value of y / y23, where 7/23 is the 
three-jet resolution parameter, using only phase-space weights. We find this average to be 
approximately equal to 0.3. Since y/y^s is, roughly, the relative transverse momentum of 

production in the range 3 < lnj/^ < 7. 
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the hardest branching in the clustering history, we select fx = 0.3-^/i as the default choice 
for the renormalization scale of a s that we use to describe the five-jet production. 

With this choice of the renormalization scale, we compare in fig. |3] ALEPH LEP1 
data for cr^do^/dlny^ and i?5 with our leading and next-to-leading results. We use 
the value of strong coupling constant a s (Mz) = 0.130 for leading order computations and 
a s (Mz) = 0.118 for next-to-leading order computations. While it is not customary to 
change the value of the strong coupling constant from one order in perturbation theory to 
the other in applications of QCD to e + e~ physics, it is done routinely in the context of 
hadron collider physics. Our choice of the leading-order value for a s is motivated by fits of 
parton distribution functions and of the strong coupling constant described in ref. [65]; our 
next-to- leading order value for a s is close to the world average [66,67]. We do not apply 
hadronization corrections at this stage. In order to assess the perturbative uncertainty, a 
scale variation by a factor of two around the default scale fi = 0.3 Mz is performed. A close 
inspection of the two plots shown in fig. || reveals that the most important effect of the NLO 
QCD corrections is the reduction in the uncertainty related to the renormalization scale de- 
pendence. The renormalization scale uncertainty is reduced from [—30%, +45%] at leading 
order to [—20%, +25%] at next-to-leading order. When leading order and next-to-leading 
order results are compared for /i = 0.3 Mz, the QCD corrections increase the leading order 
predictions by 10 — 20% 5 . The agreement between NLO QCD results and ALEPH data is 
very good for both observables considered. However, for the a^do^s/dlny^ 1 distribution 
systematic differences appear for hiy± 5 > 5.2, whereas the R§ data can be described by 
fixed-order QCD without hadronization corrections all the way up to lny~^ t = 6.5. 

4. The strong coupling constant from five-jet observables 

It follows from the discussion in the previous Section that the use of the world-average 
value of the strong coupling constant results in good agreement between parton-level NLO 
predictions, and LEP1 data for the five-jet resolution parameter and the five-jet rate. 
Therefore, we can turn this consideration around, and extract the value of the strong 
coupling constant from these two five-jet observables at LEP1 and LEP2. To combine the 
values of the strong coupling constant extracted from different observables and at different 
energies, we use the procedure advocated by the LEP QCD working group [68]. Technical 
details of the fit procedure are reported in the Appendix. 

We consider ALEPH data [16] for R§ and cruder /dy^s, as measured at LEP1 (y/s = Mz), 
and at LEP2 (^/i = 183, 189, 200, and 206 GeV). We point out that data for a^ ]da/dy 45 
are not available at ^fs = 200 GeV. Similar to other precision-QCD studies at LEP, we do 
not consider LEP2 data taken at ^/s < 183 GeV. As explained above, we use SHERPA 
for the computation of hadronization corrections. By default, hadronization in SHERPA 
is performed with the cluster model [64], which we therefore adopt as our default as well. 
Incidentally, we observe that it gives smaller hadronization corrections in the kinematic 
range of interest, than the Lund string hadronization model [63], which is also available in 

5 We note that had we used the same value of the coupling constant at LO as at NLO, as is usually done 
in e + e~ calculations, NLO corrections would have been much larger, 45 — 60%. 
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LEP1, hadr. 

°tot dcr / d 2/45, R5 


LEP1, no hadr. 

a tol da / d y^, R5 




+0.0002 


+0.0002 


stat. 






-0.0002 


-0.0002 


syst. 


+0.0027 


+0.0027 


-0.0029 


-0.0029 


pert. 


+0.0062 


+0.0068 


-0.0043 


-0.0047 


fit range 


+0.0014 


+0.0005 


-0.0014 


-0.0005 


hadr. 


+0.0012 








-0.0012 




a s (M z ) 


0.1159 +0 -° 070 
-0.0055 


0.1163 +a ° 073 
-0.0055 



Table 2: Values of the strong coupling constant a s (Mz) obtained from fits to ALEPH LEP1 data 
for cr to Jdcr/d?/45 and R 5 . NLO QCD predictions are used. Hadronization corrections are estimated 
with SHERPA. Default fit ranges are 3.8 < -lny 45 < 5.2, and 4.0 < -lnt/ cut < 5.6. See the text 
for details. 

SHERPA. We use the Lund string model to estimate systematic uncertainties related to 
hadronization effects. 

Since we use fixed-order perturbative results and do not perform any resummation, it 
is not possible to describe the data in the full kinematic ranges studied by experiments. 
This feature makes the choice of the kinematic range used in the fit an important but, 
unfortunately, somewhat a subjective issue. In general, we attempt to take the fit range 
as large as possible, with the condition that our computations are reliable and that the 
data quality is good. In the determination of the central value of a s at LEP1, we consider 
3.8 < — In 2/45 < 5.2 (7 data points) for the five-jet resolution parameter distribution, and 
4.0 < — lny cut < 5.6 (8 data points) for R$. In order to estimate the error on a s related 
to our choice of the fit range, we extract the value of ot s by performing a second fit, with 
larger ranges 3.4 < — lnj/45 < 5.6 (11 data points) for the five-jet resolution parameter, 
and 3.4 < — lny cut < 6.0 for i?5 (13 data points). The difference between the values of 
a s obtained in the two fits is called the "fit range" error; it is supposed to quantify the 
uncertainty on a s due to the choice of the data points included in the fits. 

At LEP2 the situation is different. Firstly, data are given with a coarser binning and, 
secondly, large fluctuations arc present in experimental results at small values of hiy - ^ 1 
and lny^ (for example, for some center-of-mass energies the corresponding observables 
are not even monotonic). Because of this, we decided to exclude those data points from 
our fits, effectively reducing the fit ranges. We note that those data would have had a 
modest impact on the final result anyhow, because they are affected by fairly large errors. 
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LEP2, no hadr. LEP2, no hadr. LEP2, no hadr. 





°"tot d °7 d y45 


R 5 


°tot dcr / d y45, #5 


stat. 


+0.0020 
-0.0022 


+0.0022 
-0.0025 


+0.0015 
-0.0016 


syst. 


+0.0008 
-0.0009 


+0.0012 
-0.0012 


+0.0008 
-0.0008 


pert. 


+0.0049 
-0.0034 


+0.0029 
-0.0020 


+0.0029 
-0.0020 


fit range 


+0.0038 
-0.0038 


+0.0030 
-0.0030 


+0.0028 
-0.0028 


a s {M z ) 


0.1189 +a0 ° 66 
-0.0057 


0.1120 +a ° 050 
-0.0047 


0.1155 +00044 
-0.0039 



Table 3: Values of the strong coupling constant a s (Mz) obtained from fits to ALEPH LEP2 
data with E cm > 183 GeV for cr t ~Jdtr /dj/45 and R5. NLO QCD predictions are used. Hadronization 
corrections are not included. Default fit ranges are 4.8 < — In j/45 < 6.4, and 2.1 < — log 10 y cu t < 2.9. 
See the text for details. 

We use 4.8 < — In 2/45 < 6.4 for the five-jet resolution parameter (2 data points per y^); 
and 2.1 < — log 10 y cut < 2.9 for R§ (4 data points per \fs), to find the central values of a s . 
In order to estimate the fit-range error, we employ 4.8 < — In j/45 < 5.6 (1 data point per 
Y^i), and 2.1 < — log 10 y cut < 2.5 (2 data points per -^/s), since the choice of these ranges 
leads to the largest changes in the values of the strong coupling constant compared to the 
a s values obtained from fitting with the default ranges. 

The results of our fits to LEP1 data are shown in table |2| The agreement between 
the two values of a s (Mz) extracted with and without hadronization corrections is impres- 
sive; the difference is completely negligible compared to the overall uncertainties. This 
result could have been anticipated by inspecting fig. |2[ which shows that, in the fit region, 
hadronization corrections are small, in particular when the default SHERPA choice, the 
cluster model, is used. We note that if we use the hadronization corrections as given by con- 
ventional HERWIG, PYTHIA, or ARIADNE, without matching them to high- multiplicity 
matrix elements, the picture changes drastically and the values of a s (Mz) extracted with 
or without hadronization corrections are quite different from each other. We also note 
that the overall errors of the two results given in table || are slightly smaller when in- 
cluding hadronization corrections. This is due to a marginally better description of the 
data in the central region of the fit range - which leads to smaller value of a s (Mz) and 
thus to smaller perturbative errors. However, the error reduction is partially compensated 
by the degradation of the fit quality when including larger values of In yj 5 x , In y~* t , where 
hadronization corrections increase. This feature leads to larger fit-range error compared to 
the no-hadronization case. Note also that if we extract the values of a s (Mz) by fitting 6 the 

6 Hadronization corrections are included. 
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five-jet resolution parameter distribution and R§ separately, we obtain 0.11681^^ and 
0.1151^QQQ5g respectively. These values are consistent with the result of the combined fit 
shown in table || but have slightly larger errors. From table ||, it is clear that the sensitivity 
of the five-jet observables to a s is very high, as illustrated by the tiny statistical errors. 
This sensitivity is ultimately related to the high power of a s that enters the five-jet ob- 
servables. In spite of this, the overall error is not particularly small, since the perturbative 
uncertainty is still quite sizable at this order in perturbative QCD. 

Compared to LEP1, there are important differences when we extract a s by fitting to the 
LEP2 data. Firstly, because hadronization corrections are negligible at LEP1, and because 
these corrections decrease with energy, we do not consider them for LEP2. Secondly, for the 
reasons explained above, we do not consider the data points at small values of In y^ 1 and 
lny~^ t . This fact, combined with coarser binning of data, pushes us to the region of 2/45 that 
may be affected by large logarithms of the resolution parameter. As a result, we find larger 
fit-range errors at LEP2 than at LEP1. The statistical errors are also much larger at LEP2 
than at LEP1, as one expects given the luminosities collected. On the other hand, since the 
effective strong coupling is smaller at LEP2, the perturbative uncertainty affecting five-jet 
observables decreases, making the a s extraction at LEP2 competitive with that done at 
LEP1. This effect is particularly strong for R§, for which the extracted central value of 
a s is slightly smaller than for the five-jet resolution parameter. In table || we present the 
a s values obtained by fitting separately the five-jet resolution parameter and R5 at LEP2, 
since they differ from each other by a larger amount than at LEP1. Still, both values 
are within one standard deviation from the strong coupling constant that we obtain by 
performing a simultaneous fit to the two observables. We take the latter value, given in 
the third column of table |3|, as our best determination of a s from LEP2 data. 

We obtain our final estimate of the strong coupling constant by combining the values 
of a s (Mz) extracted from LEP1 and LEP2 data. We assume that the statistical and sys- 
tematic errors of the two results are not correlated (an assumption which is strictly correct 
for the former, and a very good approximation for the latter), while the perturbative errors 
are considered to be fully correlated. The correlation of the perturbative uncertainties is 
due to the fact that we estimated them by varying the renormalization scale, which results 
in changes of the cross sections whose pattern is independent of the center-of-mass energy. 
It is quite likely that a more sophisticated approach to estimating perturbative errors (see 
e.g. ref. [69]) will result in a smaller uncertainty on a s . Hence, the procedure that we 
employ in this paper is rather conservative. Using the results of tables Q and p|, we finally 
obtain 

a s (M z ) = 0.11561^034 ■ (4-1) 

We note that if we perform the fit to both LEP1 and LEP2 data simultaneously, we obtain 
a s (M z ) = 0.1156lo:oo4l' in perfect agreement with eq. ([O]). 

The value of a s (Mz) that we extract from five-jet observables at LEP can be compared 
with other recent determinations of this quantity, shown in table ||. We see that both the 
central value of a s and its error, obtained from fitting five-jet observables, compare well 
with other determinations. On the other hand, it is interesting that as(Mz) in eq. ([O]) 
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Observable 


aJMz) 


Ref. 


r decays 


0.1197 ± 0.0016 


[67] 


X decays 


n i in -i- n nn^ 


r7ni 


3 jet observables 


0.1224 ±0.0039 


[44] 


jets in DIS 


0.1198 ±0.0032 


[71] 


DIS 


0.1142 ± 0.0021 


[72] 


thrust 


0.1135 ±0.0011 


[73] 


lattice 


0.1183 ±0.0008 


[74] 


EW fits 


0.1193 ±0.0028 


[75] 


world average 


0.1184 ±0.0007 


[67] 


e + e~ — > five jets 


0.1156 ±0.0038 


this paper 



Table 4: Summary of selected determinations of the strong coupling constant a s (Mz). We have 
averaged the ± errors shown eq. (4.1) for the determination of a s reported in this paper. Not all 
results shown in this table are included in the world average. 



is lower than the world average. It is peculiar that a number of recent determinations of 
a s arrived at a similar conclusion. 



5. Conclusions 

In this paper we study the production of five jets in e + e~ annihilation at LEP1 and LEP2. 
We improve the perturbative QCD predictions for five-jet observables, l/o"totd<r/dy45 and 
i?5, by computing the NLO QCD corrections. For suitably chosen renormalization scales, 
such corrections are between ten and twenty percent 7 . They reduce the scale uncertainty 
by about a factor of two with respect to the LO predictions, and lead to a better agreement 
between theoretical predictions and experimental data. 

We point out that hadronization corrections computed with event generators whose 
showers are not matched to high-multiplicity matrix elements (such as out-of-the-box HER- 
WIG, PYTHIA, and ARIADNE) are large and uncertain. For this reason, we believe it is 
important to describe five-jet observables in a way that incorporates high-multiplicity tree- 
level matrix elements. This is provided by the event generator SHERPA, which implements 
the CKKW procedure for matching tree- level matrix elements to parton showers. In this 
way, an improved description of five hard, well-separated partons is obtained, which in turn 
results in fairly small hadronization corrections in the range where fixed-order perturbative 
results are most reliable. 

We extract the strong coupling constant from the distributions of the five-jet reso- 
lution parameter and the five-jet rate, as measured at LEP1 and LEP2 by the ALEPH 

7 Note that this statement is only valid if the leading order result is calculated with ces(Mz) = 0.130, 
which is much larger than the value of the strong coupling constant given in eq. ( |4.l| ) or the world average. 
Changing the value of the strong coupling constant from one order in perturbation theory to the other is 
not customary in e + e~ collider physics, while it is an accepted practice in hadron collider physics. 
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collaboration. We find a s (Mz) = 0.1156^q'qq 34 , which compares well with other recent de- 
terminations of the strong coupling constant, and is somewhat lower than the current world 
average value. We stress that our treatment of the uncertainties on a s is conservative. A 
detailed knowledge of the experimental systematics, and a more sophisticated approach to 
theoretical errors will very likely lead to a higher precision in the determination of a s from 
five-jet observables at LEP. 
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Appendix: details of the fit 

In this Appendix, the details of the fitting procedure are described. In the fit, we consider 
two observables - the five-jet resolution parameter distribution and the five-jet rate. These 
observables are measured at several energies; the results are available in the form of binned 
distributions. In principle, we can extract the value of the strong coupling constant from 
any of the bins but, clearly, the availability of many bins helps in decreasing the errors. The 
problem is that both the experimental and theoretical errors affecting different bins may be 
correlated, and it is important to take these correlations correctly into account, to obtain a 
proper estimate of the overall uncertainty in the determination of a s . Our procedure follows 
closely the approach of the LEP QCD working group described in ref. [68]. A possible way 
to treat theoretical uncertainties is discussed in ref. [69] , but we follow a simpler approach 
that is described below. 

The set of all available bins, for a chosen range of 2/45 and y C ut, is a set of observables 
La from which values of the strong coupling constant can be determined. Each member 
of Lq is represented by three numbers Oi = [Xj, af tat , <r| yst ], where X i} and a ^ tat ' s y st a re 
the central value and the statistical and systematic uncertainties respectively, for a given 
observable in the bin i. If uncertainties are asymmetric, Oi is a collection of five numbers, 
and the discussion below applies to positive and negative errors separately. In what follows, 
we use a s ± 5a s as the shorthand notation instead of the full one a s + ft as . 

As we already mentioned, each of the observables from the list Lq can be used to 
determine the value of the strong coupling constant. This is done by solving the equation 



Ti Hi — Ei, 



(5.1) 



where Tj is the (parton level) theoretical prediction, Hi is the hadronization correction, 
and Ei is the experimental value for the bin i. The theoretical prediction Tj depends on 
a s (Mz) and the renormalization scale \x. As discussed in the text, for the hadronization 
correction we can use either the cluster or the Lund string model and we choose the former 
as our default. We write the value of the strong coupling constant at Mz, obtained by 
solving eq. ( |5.1| ), as 

a* = al ± <5c4' stat ± 5al' syst ± Jc4' scale ± 5af iadv , (5.2) 

where a l s is the central value. The central value and the errors in eq. ( |5.2| ) are obtained in 
the following way: 

1) the central value o| is obtained by solving eq. (pJj) for a s with /i = [io, no = 0.3y/s, 
Ei = Xi (i.e., the central data value in the relevant bin), and the cluster model for 
hadronization. The results without hadronization corrections are obtained by simply 
setting Hi = 1. 

2) 5a l s ' stat and 5a l s ' syst are obtained by solving eq. ( |5.1| ) for a s with the same settings 
as in item 1), except that Ei = Xi ± cr l ' stat or Ei = Xi ± a l ' syst . The differences 
between the values of a s obtained in this way, and the central value ot s , are ±5as Stat 
and ±5al' syst . 

3) 5af cale is obtained by solving eq. ( fP| ) for a s with the same settings as in item 1), 
except that \x = 0.5 no and \i = 2^0 are used. The differences between the values of 
a s obtained in this way, and the central value a*, are ±(5as' scae . 

4) <5oi' hadr is obtained by solving eq. (|5.1| ) for a s with the same settings as in item 1), 
except that the Lund string model is used for hadronization. We define ±<5as' hadr = 
±|a s — al \. We note that this is a conservative choice, since clearly a s — a* is either 
positive or negative. 

The result of this procedure is a set of values of the strong coupling constants, a\, with the 
corresponding errors. They need to be combined to obtain the average value of a s . To this 
end, it is necessary to construct the covariance matrix, which we define as the sum of the 
covariance matrices for statistical, systematic, perturbative, and hadronization errors. If 
we denote generically by 5a l s one of these four errors, the corresponding covariance matrix 
is 

Vij = dij {5a l s f + (1 - Sij) dj5a l M , (5.3) 

where Cjj is the statistical correlation between a\ and oP s (see later). The covariance matrix 
is used to calculate the average value of the strong coupling constant and its error in a 
standard way. We compute the weights 

N N 

w i = ^(V-%/ (V-'hi, l<i<N, N = dim(V), (5.4) 

j=l k,l=l 



- 18 - 



and obtain the estimate of the average of the strong coupling constant and of its error 



N 



a s = '^2w i a l s , (5.5) 



i=l 

N 



a 2 (a s ) = WiVijWj . (5.6) 

It is easy to see that if two or more errors that enter the definition of the covariance matrix 
are close numerically, and the absolute values of the corresponding statistical correlations 
Cij is close to one, eq. (|5.3| ) may lead to pathological results, since the weights Wi tend 
to grow large in absolute value and to have opposite signs. To avoid this, the LEP QCD 
working group [68] adopts the formula 

Vij = Sij (6a\f + (1 - dij) min{(K) 2 , {5cP s ) 2 } , (5.7) 

which is essentially equivalent to taking the largest possible Cij in eq. (^) , that still leads 
to non-pathological weights W{. Clearly, eq. ( p.7[ ) is an overestimate of the correlation if 
the actual Cij is small. For this reason, in our fit we use a slightly modified formula for 
the covariance matrix 

Vij = 5ij (5al) 2 + (1 - S^) min {(5al) 2 , (5ai) 2 , Cyfajtfoj} . (5.8) 



As we explain below, for the observables and errors that we consider, eq. (|5.8|) coincides 
with eq. (5.7) in all cases, except for the statistical correlation between a^da /dy^ and 



In summary, we take the covariance matrix to be 



V = V stat + V syst + y scale -f y hadr (5 9) 



where each of the terms on the right-hand side of eq. ( |5.9D is constructed according to 

A r |\ ■ r i.stat r i,syst r i, scale j r i,hadr , ■ i • r j.i xr j* 1 

eq. ( |5,a| ) using oa s , oa s , oa s , and oa s respectively. As far as the off-diagonai 



terms of the various V matrices are concerned, we have assumed what follows: 

• Statistical errors are uncorrelated between different center-of-mass energies. At a 
given center-of-mass energy, <r t ~ \ da/dy^ data are uncorrelated, data are fully 
correlated, and R$ data are correlated with a^ \da / dy^ ones for y cnt < 2/45. We have 
explicitly computed the coefficient Cij relevant to the cr t ~Jd<7 / dy^ — R5 correlation, 
and have used them in eq. ( |5.8[) . 



Since we do not have detailed information about correlations of systematic uncertain- 
ties in ALEPH data, we assume conservatively that all systematic errors at a given 
center-of-mass energy are fully correlated. We also assume that systematic errors are 
completely uncorrelated between LEP1 and LEP2, but that they are fully correlated 
in the measurements performed at different LEP2 energies. 
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• Perturbative errors are assumed to be fully correlated, for all observables and energies. 
See also the main text for a comment on this point. 

• Hadronization errors, that we compute only at LEP1, are assumed to be fully corre- 
lated. 

Finally we point out that in the computation of the central value of a s , according to 



eq. (5.5), we neglect the off-diagonal entries of y scale and y hadr [68]. On the other hand, we 



include all the off-diagonal entries in the computation of the standard deviation according 



to eq. (5.6). In fact, in the presence of errors numerically very close to each other (which is 
the case for the perturbative and hadronization errors) , the result for the average value of a s 
tends to assume the value of the input with the smallest error, which is again an artifact of 
the combination procedure. We have checked that the value of the strong coupling constant 
that we obtain in this way is statistically fully compatible with the result we would have 
obtained by considering all correlations when determining the central value of a s . 
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